Rozważmy dla ustalenia uwagi nasz problem adwekcji-dyfuzji (obliczenia propagacji zanieczyszczeń na danym obszarze za pomocą dyfuzji oraz wiatru)
\( b(u,v)=l(v) \\ b(u,v) = \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y)}{\partial x} dxdy \int_{\Omega } \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ \int_{\Omega} \epsilon \frac{\partial u(x,y)}{\partial x } \frac{\partial v(x,y)}{\partial x } dxdy +\int_{\Omega} \epsilon \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \\ l(v) = \int_{\Omega} f(x,y) v(x,y) dxdy \)
Oznaczamy przez \( U \) zbiór funkcji w którym zawiera się rozwiązanie problemu adwekcji-dyfuzji \( u \in U \).
Następnie przez \( V \) oznaczmy zbiór funkcji które używamy do testowania w naszym sformułowaniu wariacyjnym (zwanym też słabym).
W rozdziale tym opiszemy za pomocą trochę abstrakcyjnego języka matematycznego jak wyprowadza się metodę minimalizacji reziduum [1].
Nasz oryginalny problem zapisać można tak: szukamy funkcji \( u \in U \) takiej że \( b(u,v)=l(v) \) dla wszystkich funkcji testujących \( v\in V \).
Matematycznie można oznaczyć lewą stronę \( Bu \) a prawą stronę \( l \) wprowadzając pewne operatory matematyczne (matematycy powiedzą że \( B : U \rightarrow V' \) oraz \( \langle Bu, v \rangle = b(u,v) \).
Wówczas nasz problem sprowadza się do znalezienia pola skalarnego \( u \) takiego że
\( Bu-l=0 \). Wyrażenie \( r=Bu-l \) nazywa się reziduum.
Jeśli mamy jedynie rozwiązanie przybliżone, na przykład oznaczone \( u_h \)
wówczas nasze reziduum będzie wyrażało błąd, im większe będzie reziduum tym błąd rozwiązania przybliżonego będzie większy \( r_h=Bu_h-l \). W jaki sposób zmierzyć błąd rozwiązania, czyli wielkość reziduum?
Żeby zmierzyć jak duże jest reziduum musimy wprowadzić jakąś normę (czyli miarę wielkości) \( r_h = \| Bu_h-l \| \). Mogę teraz powiedzieć że szukamy takiego rozwiązania naszego problemu \( u_h \) żeby reziduum (czyli błąd) był jak najmniejszy. Matematycznie zapisuje się ten warunek w postaci minimalizacji normy
\( u_h = \displaystyle{argmin}_{w_h \in U } {1 \over 2} \| Bw_h - l \|_{V' }^2 \)
Przed normą dopisuje się \( 1\over 2 \) oraz podnosi się normę do kwadratu.
Problem praktyczny jest taki że nie wiemy jak policzyć normę z różnicy \( \|Bw_h-l\| \). Norma ta mierzy bowiem odległości w abstrakcyjnej przestrzeni matematycznej \( V' \) (matematycy powiedzą że przestrzeń ta jest przestrzenią dualną do przestrzeni funkcji testujących).
Żeby rozwiązać ten problem, wykonujemy trick matematyczny polegający na rzutowaniu normy z przestrzeni \( V' \) do przestrzeni testującej \( V \). Wprowadza się operator przerzucający przestrzeń \( V \) w przestrzeń \( V' \). Operator ten nazywamy operatorem Riesza \( R_V \colon V \ni v \rightarrow (v,.) \in V' \). Operator odwrotny natomiast przerzuca przestrzeń \( V' \) z powrotem w przestrzeń \( V \), czyli \( R_V^{-1} \colon V' \ni (v,.) \rightarrow v \in V \). Tak więc jeśli użyjemy odwrotnego operatora Riesza, przerzucimy nasz problem do przestrzeni funkcji testujących \( V \). Mamy więc
\( u_h = \displaystyle{argmin}_{w_h \in U_h} {1 \over 2} \| {\color{red}{R_V^{-1}}} (Bw_h - l) \|_{\color{red}{V} }^2 \). Szukamy więc punktu w którym funkcja osiąga minimum. Taki punkt osiągany jest wtedy kiedy pochodne tej funkcji są równe zero. Ponieważ mamy tutaj funkcjonał, matematycy powiedzą że szukamy takiego \( u_h \) w którym pochodne Gateaux są równe zero we wszystkich kierunkach
\( \langle R_V^{-1} (Bu_h - l), R_V^{-1}(B\, w_h) \rangle_V = 0 \quad \forall \, w_h \in U_h \) gdzie \( w_h \) oznacza wszystkie możliwe "kierunki" ze skończenie wymiarowej przestrzeni funkcji \( U_h \) (w której przy okazji szukamy rozwiązania naszego problemu adwekcji-dyfuzji).
Formalnie rzecz ujmując, w tym miejscu musimy skonstruować skończenie wymiarową przestrzeń funkcji \( U_h \) w której będziemy szukać naszego rozwiązania. Robimy to oczywiście w naszym podręczniku stosując funkcje bazowe B-spline, rozpinając je zgodnie z wektorami węzłów wzdłuż osi x oraz wzdłuż osi y.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4].
Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline \( \{ B_{i,2}(x) \}_{i=1,...,6 } \) oraz \( \{B_{j,2}(y)\}_{j=1,...,6 } \). Następnie utworzymy z nich iloczyny tensorowe \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \).
Zauważmy że nasz obszar \( \Omega \) rozpięty jest na kwadracie \( [0,1]^2 \), podobnie jak nasze 6*6=36 funkcji bazowych co wynika z definicji wektora węzłów [0 0 0 1 2 2 3 4 4 4].
Nasza przestrzeń aproksymacyjna jest więc rozpięta przez iloczyny tensorowe naszych jednowymiarowych funkcji B-spline
\( U_h = gen \{ B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \} \).
Jest więc ona przestrzenią 6*6=36 wymiarową, Mamy więc 36 kierunków (wielomianów bazowych B-spline) w których liczyć będziemy pochodne kierunkowe i będziemy je przyrównywać do zera.
Oznaczamy
\( r=R_V^{-1}(Bu_h-l) \)
i wówczas nasz problem minimalizacji reziduum zapisać można w postaci
\( \langle r , R_V^{-1} (B\, w_h ) \rangle = 0 \quad \forall \, w_h \in U_h \)
co jest równoznaczne
\( \langle Bw_h, r \rangle = 0 \quad \quad \forall w_h \in U_h. \)
Dorzucamy teraz drugie równanie, korzystając z definicji reziduum \( r=Bu_h-l \), przemnażam definicje reziduum przez funkcje testujące z przestrzeni funkcji testujących i dostaje równanie:
\( (r,v)_V=\langle B u_h-l,v \rangle \quad \forall v\in V. \)
Mam więc układ równań: znajdź reziduum w nieskończenie wymiarowej przestrzeni funkcji testujących \( r\in V \) oraz rozwiązania w skończenie wymiarowej przestrzeni funkcji \( u_h \in U_h \)
\( (r,v)_V - \langle Bu_h-l ,v \rangle = 0 \quad \forall v\in V \\ \langle Bw_h,r\rangle = 0 \quad \forall w_h \in U_h \)
Gdybyśmy byli wstanie rozwiązać to równanie, to dostali byśmy najlepszą możliwą aproksymacje naszego problemu w przestrzei funkcji aproksymujących \( U_h \). Niestety nie jest to możliwe, ponieważ nieskończona przestrzeń funkcji testujących generuje nam nieskończenie wiele równań które należało by rozwiązać. Musimy więc wybrać drugą inną bazę do przybliżenia przestrzeni testowej. Tutaj pojawia się pierwsza istotna różnica metody minimalizacji reziduum z tradycyjną metodą elementów skończonych. Otóż w metodzie minimalizacji reziduum mamy dwie różne przestrzenie, jedną do aproksymacji rozwiązania, drugą do testowania.Oczywiście przestrzeń testującą również rozpinamy za pomocą funkcji B-spline. Na przykład możemy użyć takiego samego patcha elementów jak dla przestrzeni aproksymującej (ale nie jest to wymagane!), albo podnieść stopień funkcji B-spline.
Precyzujemy jak wyglądać będzie nasz patch elementów oraz nasze funkcje bazowe B-spline rozpięte na elementach, podając wektor węzłów wzdluż osi \( x \) oraz wektor węzłów wzdłuż osi \( y \). Odsyłamy tutaj do rozdziału trzeciego w którym opisany jest sposób definiowania funkcji bazowych za pomocą wektorów węzłów i wielomianów B-spline.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 0 11 2 2 3 3 4 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 0 1 1 2 2 3 3 4 4 4 4]. Generalnie przestrzenie testujące muszą być większe niż przestrzeń aproksymacyjna. Możemy uzyskać ten efekt zwiększając liczbę elementów lub podnosząc stopień funkcji B-spline, lub redukując ciągłość bazy funkcji B-spline, lub mieszając te trzy podejścia razem.
Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline \( \{ \tilde{B}_{k,3}(x) \}_{k=1,...,10 } \) oraz \( \{\tilde{B}_{l,3}(y)\}_{l=1,...,10 } \). Następnie utworzymy z nich iloczyny tensorowe \( \tilde{B}_{k,l;3}(x,y)=\tilde{B}_{k,3}(x)B_{l,3}(y),k,l=1,...,10 \).
Nasza przestrzeń testujące to
\( V_h = gen \{ \tilde{B}_{k,l;3}(x,y)=\tilde{B}_{k,3}(x)B_{l,3}(y),i,j=1,...,10 \} \).
Mając dyskretną przestrzeń testującą, uzyskujemy wreszcie nasz dyskretny (czyli skończenie wymiarowy) układ równań który będziemy chcieli rozwiązać.
Szukamy \( ( r_m, u_h )_{ V_m \times U_h } \)
\( (r_m,v_m)_{V_m} - \langle B u_h-l,v_m \rangle = 0 \quad \forall v\in V_m \\\langle Bw_h,r_m\rangle = 0 \quad \forall w_h \in U_h \)
Pojawia się tutaj tak zwany iloczyn skalarny (wewnętrzny) w dyskretnej przestrzeni testującej \( (*,*)_{V_m} \).
Żeby odgadnąć jaki iloczyn skalarny użyć w naszym problemie, musimy popatrzeć na nasze oryginalne równanie zapisane w formie słabej
\( b(u,v)=l(v) \\ b(u,v) = \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y)}{\partial x} v(x,y) dxdy+ \int_{\Omega } \beta_y(x,y) \frac{\partial u(x,y)}{\partial y }v(x,y) dxdy +\\ \int_{\Omega}\epsilon \frac{\partial u(x,y)}{\partial x } \frac{\partial v(x,y)}{\partial x } dxdy +\int_{\Omega} \epsilon \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \\ l(v) = \int_{\Omega} f(x,y) v(x,y) dxdy \)
i zobaczyć jakiej formy są funkcje testujące \( v \) w formie \( b(u,v) \). Widzimy że mamy tutaj doczynienia z pochodnymi cząstkowymi \( \frac{\partial v(x,y)}{\partial x } \), \( \frac{\partial v(x,y)}{\partial y } \) oraz z samą funkcją \( v(x,y) \). Nasz produkt wewnętrzny powinien więc zawierać pochodne i wartości funkcji
\( (u,v)_{V_m} = \int_{\Omega} u(x,y) v(x,y) dxdy + \int_{\Omega} \frac{\partial u(x,y)}{\partial x } \frac{\partial u(y,y) }{\partial x } dxdy + \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial u(y,y)}{\partial y } dxdy \)
Podsumowując, w metodzie minimalizacji reziduum, musimy zdefiniować osobną przestrzeń aproksymacyjną, oraz osobną (większą) przestrzeń testującą, oraz wybrać iloczyn wewnętrzny przestrzeni testującej. Dostajemy wówczas układ równań w którym niewiadome to nasze rozwiązanie \( u \) oraz dodatkowo reziduum \( r \).
To jak dobrze działać będzie metoda minimalizacja reziduum zależy w dużej mierze od naszego wyboru przestrzeni testującej oraz iloczynu wewnętrznego. Jeśli iloczyn wewnętrzny będzie dostatecznie dokładny, oraz przestrzeń testująca będzie dostatecznie duża, wszystko będzie idealnie działać, i dostaniemy najlepsze możliwe do uzyskania rozwiązanie w przestrzeni aproksymacyjnej \( U_h \) (ale nie lepsze niż pozwala na to przestrzeń aproksymacyjna). Matematyczne uzasadnienie jest takie, że problem minimalizacji reziduum z nieskończenie wymiarową przestrzenią testową realizuje warunek stabilności inf-sup ze stałą równą 1. Jeśli przestrzeń testującą zawęzimy do przestrzeni skończenie wymiarowej, to warunek inf-sup może nie być już idealnie realizowany w tej przestrzeni. Zwiększanie przestrzeni testującej będzie nam poprawiać realizowalność warunku stabilności inf-sup.
Nasz układ równań do rozwiązania w metodzie minimalizacji reziduum wyglada następująco
\( \left( \begin{matrix} G & B \\ B^T & 0 \\ \end{matrix} \right) \left( \begin{matrix} r \\ u \end{matrix} \right)= \left( \begin{matrix} f \\ 0 \end{matrix} \right) \)
Zobaczmy jak wygląda nasz układ równań na przykładzie konkretnych przestrzeni aproksymacyjnej i testującej.